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ABSTRACT 

Context. In the recent years, MRI-driven turbulent transport has been found to depend in a significant way on fluid viscosity v and 
resistivity rj through the magnetic Prandtl number Pm = v/rj. In particular, the transport decreases with decreasing Pm; if persistent 
at very large Reynolds numbers, this trend may lead to question the role of MRI-turbulence in YSO disks, whose Prandtl number is 
usually very small. 

Aims. In this context, the principle objective of the present investigation is to characterize in a refined way the role of dissipation. 
Another objective is to characterize the effect of linear (channel modes) and quasi-linear (parasitic modes) physics in the behavior of 
the transport. 

Methods. These objectives are addressed with the help of a number of incompressible numerical simulations. The horizontal extent 
of the box size has been increased in order to capture all relevant (fastest growing) linear and secondary parasitic unstable modes. 
Results. The major results are the following: 

i- The increased accuracy in the computation of transport averages shows that the dependence of transport on physical dissipation 
exhibits two different regimes: for Pm < 1, the transport has a power-law dependence on the magnetic Reynolds number rather than 
on the Prandtl number; for Pm > 1, the data are consistent with a primary dependence on Pm for large enough (~ 10 3 ) Reynolds 
numbers. 

ii- The transport-dissipation correlation is not clearly or simply related to variations of the linear modes growth rates. 

iii- The existence of the transport-dissipation correlation depends neither on the number of linear modes captured in the simulations, 
nor on the effect of the parasitic modes on the saturation of the linear modes growth. 

iv- The transport is usually not dominated by axisymmetric (channel) modes. 
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1. Introduction 

Disks evolve on time-scales that are orders of magnitudes 
smaller than expected from microphysical transport processes, 
and various suggestions have been made over the years to ex- 
plain this discrepancy. Turbulent transport, in particular, has fig- 
ured among the leading candidates since the inception of the a- 
disk paradigm, and a number of hydrodynamic and MHD turbu- 
lent transport mechanisms have been proposed in the literature. 

On the hydrodyn amic side, subcritical turbulence 
( Richard & Zahn 1999 and references therein), if present, 
is app arently too inefficient dLesur & Longarettill2005t1ji et al.l 
2006). Convection was up to now found too inefficient and to 
transport angular momentum in the wrong direction (Cabot 
1996; Stone & Balbusl Il996h . but a recent reinvestigation of 
the problem indicates that this might be an artifact of these 
simulations being perfor med too close to the stability threshold 
dLesur & Ogilvi e 2010). Two-dimensional weak turbulence 
driven by small-scale, incoh erent gravitation al instabilities 
(density waves) is an option (Gammie 1996). Alterna tively, 
the baroclinic instability (Klahr & Bodenheimer 2003) may 
generate vorticity, and transport through the coupling with 
density wa ves, but its conditions of e xistence are st i ll con - 
troversial dJohnson & Gammid l2006t iPetersen et all l2007l) . 



although iLesur & Papaloizoul (1201 Oh have probably identified 
the root of this debate by pointing out the nonlinear nature of 
the instabilit y; also the resulting vortice s would be subject to 3D 
instabilities (Lesur & Papaloizou 2009). 

Balb us & Hawlevl (1 1 99 lb have proposed that the magne- 
torotational instability (MRI) is a potentially efficient source 
of turbulent transport in the nonlinear regime, an expectation 
soon borne out in numerical simulations. This instability pro- 
vides by now the most extensiv ely studied transpo rt mecha- 
nism, through loca l unstratified (Ha wlev et al.lll995l) . stratified 
(IStone et al.ll 19961) . and global dHawlevI 120001) 3D disk simu- 
lations. These initial simulations as well as the numerous ones 
following them have shown that MRI turbulence is an efficient 
way to transport angular momentum, in the presence or ab- 
sence of a mean vertical or toroidal field, with an overall trans- 
port efficiency depending on the field configuration and strength. 
However, the significant role played by microphysical dissipa- 
tion in the resolutions accessible to date had largely been u nder- 
estimated dLesur & Longarettill2007l: iFromang et al.ll2007l) . 

By now, both the field strength and dissipation dependence 
of the simulated turbulent transport have been studied to some 
extent (and only in unstratified local shearing box settings for 
the latter one). The dependence of the Shakura-Sunyaev a pa- 
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ramete r has been characterized very early on by Hawley et al. 
(1 19951) who showed that momentum transport oc pr l l 2 both for 
a net vertical or toroidal fielcQ (albeit with very different effi- 
ciencies in the two configurations), a scaling further con firmed 
in later simulations, as summarized in Pessah et al] d2007l) . 

Until recently, the effect of physical viscosity (v) and resis- 
tivity (77) on the transport had been neglected, under the implicit 
assumption that these should not matter too much once iner- 
tial turbulent scales are re solved in the simulations. However, 
iLesur & Longaretti (2007]) have shown that, in the presence of 
a mean vertical field, the MRI-driven turbulent transport did ex- 
hibit a substantial dependence on the magnetic Prandtl number 
Pm = v/77, with no clear tre nds with respect to either viscosity of 
resistivity alonfl Recently, Simon & Hawley (2009) found sim- 
ilar results in shearing boxes with a mean toroidal field instead 
of a mean vertical one. 

When the mean magnetic flux vanishes, the tra nsport behav- 
ior is more complex. The initial investigation by lHawlev et alj 
(1996J) co ncluded that the transpor t was converging to a finite 
value, but iGardiner & Stond d2005h found that the transport ef- 
ficiency was dependent on the simulation resolution. More re- 
cently, the role of the magnetic Prandtl nu mber Pm has been 
identified in this setting (Fromang et al. 2007): turbulence exists 
only for magnetic Prandtl numbers larger than about 2, which 
requires the explicit inclusion of viscous and resistive terms in 
the fluid equations for numerical simulations to correctly capture 
the physics of the problem. The disappearance of turbulence at 
low Pm, as well as the need of large enough amplitudes in the 
initial conditions at Pm > 2, indicate that the zero net flux mag- 
netized shea ring box is a subcritical sy stem rather than a linearly 
unstable one dLesur & Ogilvidl2008afbh . 

Thus it appears that in all configurations explored to date, 
the magnetic Prandtl number plays a significant role on the ex- 
istence and/or efficiency of the turbulent transport, at least at the 
resolutions accessible on present day computers (or equivalently, 
the accessible Reynolds and magnetic Reynolds numbers). This 
raises a number of issues. 

For one, the exact role played by channel modes and para- 
sitic modes is unclear. Although they exist only when a mean 
vertical field is present, they are simpler to analyze and their be- 
havior may provide insight into the generic mechanism respon- 
sible for saturation of the linear instability. Channel modes are 
the axisymmetric unstable m odes of the MRI (Balbus & Hawlevl 
ll99lHPessah & Chanll2008l) . and are often observed both in 2D 
and 3D simulations with a mean vertical field; their name derives 
from their vertically layered characteristic channel-like radial 
flow. They were quic kly recognized to be a lso nonlinear solu- 
tion of the problem bv lGoodman & XvJ d!994l) ; the same authors 
found them to be unstable with respect to a secondary instability 
(parasitic modes). A few recent papers have focused on the pos- 
sibility that the saturation of the channel mode by this parasitic 
instability might be the mechanism explaining the magnitude of 
the turbul ent transport in MRI simulations, with diver ging con- 
clusions ([PiiiihAGoodmaBllOOl; EatteretaD|200§- 

In relation to this, the role of the aspect ratio of the simula- 
tions has probably been underestimated in the past. Boxes with 
an aspect ratio R : Z = 1 : 1 do no t allow for the fastest para- 
sitic modes to grow, and Bo do et alJ (Hp08) pointed out that nar- 



1 The definition of p (ratio of gas to magnetic field pressure) is based 
on the mean field, a conserved quantity in the settings used in these 
simulations. 

2 The linear stability dependence on viscosity and resistivity is totally 
different. 



row boxes tend to overemphasize the role of the channel modes 
with respect to more horizontally extended boxes. This calls for 
a reassessment of the Prandtl number dependence of MRI-driven 
transport in horizontally extended simulation boxes with a mean 
vertical field. 

More generally, it is still unclear whether this dependence 
of the transport on physical dissipation is a consequence of the 
limited Reynolds numbers that can be achieved on present day 
computers. In particular, none of the published simulations has 
been able to capture the existence of a significant inertial range in 
the kinetic or magnetic energy spectrum, which makes it difficult 
to address this issue. The question here revolves mostly around 
the direction and locality of transfers and fluxes in Fourier space, 
and will be addressed elsewhere. 

For the time being, we focus the potential role of the channel 
and parasitic modes in the efficiency of turbulent transport. This 
is explored by numerical simulations in the shearing sheet limit, 
with a net vertical magnetic flux, and with horizontally extended 
simulation boxes. The paper is organized in the following way. 
Our numerical method, setup, and run parameters are described 
in section 2. Relevant aspects of the theory of channel modes 
are summarized in section 3. Section 4 is the core of this paper, 
and discusses our numerical results; the issues bearing on the 
resolved linear and secondary modes are also discussed there. 
The implications of these results are presented in the final section 
along with some possible future lines of work. 

2. Numerical model 

2.1. Shearing box model and equations: 

Following the initial investig ation of 3D MRI turbulent prop- 
erties by lHawlev et al. (1995), we base our simulations on the 
shearing sheet local approximation and the related shearing box 
model. Most if not all local studies of disk turbulence have been 
performed in this framework. Local simulations are unescapable 
to examine in any detail the structure and transport properties 
of MRI turbulence; indeed, even within a local model, present 
day computers are still too limited to reach the resolutions re- 
quired to understand the magnetic Prandtl issue summarized in 
the introduction, and it is certainly hopeless to tackle this prob- 
lem directly in global simulations. 

To some extent, a shearing box biases the role of the channel 
mode in turbulent transport, e.g. through the correlations intro- 
duced by the periodic boundary conditions. Note however, that, 
in purely hydrodynamic turbulence, the shearing box seems to 
capture some of the correct physical properties of actual exper- 
imental systems, such as the transition Reynolds number to tur - 
bulence as a function of rotation (see. lLesur & Longaretti 2005). 
In any case, it is very difficult to formulate a well-posed local 
problem that does no t rely on the shearing box framework . The 
reader is referred to |Hawlev et ail d 19951) . iBalbusi (l2003h and 
iRegev & Umurhanl (|2008l) for more detailed discussions of the 
properties and limitations of this model. 

MHD turbulence in discs is essentially subsonic, and we 
will work in the incompressible approximation, which al- 
lows us to eliminate local and transient density fluctuations 
as well as sound waves and density waves from the prob- 
lem. Density waves are excited in shearing box turbulence 
dHeinemann & Papaloizoul l2008b), but they appear to have lit- 
tle im pact on the turbulent transport ( Hein emann & Papaloizoul 
2008a). This has been confirmed by direct comp arison between 
incom pressible and compressible simulations dFromang et al.l 
120071) . As a consequence, we feel reasonably justified to as- 
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sume incompressibility. Explicit molecular viscosity and resis- 
tivity are included. 

The shearing box equations are well-known; we reproduce 
them here to introduce our notations. We chose a Cartesian box 
centered at r — Rq, rotating with the disc at angular velocity 
Q = Q(/?o) and having dimensions (L x ,L y ,L z ) with L, <K Rq. 
By convention here, Rq<[> = x and r - Rq — —y, leading to the 
following form of the shearing sheet equations: 



d t U + V • (U ® U) 

d,B 
V U 
VB 



- vn + V • (-B ® B) 

-2SI X U + 2D.Sye y + vAU, (1) 

V X (U X B) + r]AB, (2) 

0, (3) 

0. (4) 



where the magnetic field is measured in units of Alfven speed. 
The mean shear S = -r<9,-Q is set to a Keplerian flow value 
S = (3/2)£2. The generalized pressure term II includes both the 
gas pressure term P/po and the magnetic one B 2 /2pi). This gen- 
eralized pressure II is fixed by the incompressibility condition 
Eq. O), and computed by solving a Poisson equation. The mag- 
netic field is expressed in Alfven-speed units, for simplicity. 

The steady-state solution to these equations is the local 
Keplerian profile Uq — S ye x . Our code computes the (turbulent) 
deviations from this Keplerian profile. Defining v = U—Uq, and 
b = B - Bo, one obtains the following equations for v and b: 



d,v = -v-Vv- Vn + B Vb- Syd x v 
+(2Q - S)v y e x - 2Q.v x e y + vAv, 

d,b = -v-Vb + B- Vv 

-Syd x b + Sb y e x + rjAb, 

V • v = 0, 

V • b = 0. 



(5) 

(6) 
(7) 
(8) 



The boundary conditions associated with this system are pe- 
riodic i n the x and z directio n and shearing-periodic in the y di- 
rection (lHawlev et al.ll 19951) (fort; and 6). 

Following Hawley et al.l d 1995b . one can integrate the induc- 
tion equation (|6]l over the volume of the box, leading to: 



d(B) 

at 



S(B y )e x , 



(9) 



where () denotes a volume average. Therefore, the mean mag- 
netic field is conserved, provided that no mean radial field is 
present. In this work, a mean vertical field Bq is imposed, and 
conserved by virtue of Eq. ((9]). 

The numerical resolution makes use of a spectral Galerkin 
representation of equation s ©-(H) in the sheared frame (see 
iLesur & Longa rettil 120051) . In this frame, the shearing-sheet 
boundary conditions are transformed into perfectly periodic 
boundary conditions, and Fourier transforms can be used in all 
three directions. Moreover, this decomposition allows us to con- 
serve magnetic flux to machine precision (the total magnetic flux 
created during one simulation is typically 10~ n ). Equations (0 
and <[8J ar e enforced to machine precision using a spectral pro- 
jection dPevrejl2002h . The nonlinear terms are computed with a 
pseudospectral method, and aliasing is prevented using the 3/2 
rule. 

The sheared frame representation of the spectral domain 
would eventually produce a mismatch between the computed 



Fourier components and the physically relevant ones. To circum- 
vent this p roblem, we use a remap m ethod similar to the one de- 
scribed bv lUmurhan & Regevl d2004h . This routine redefines the 
sheared frame every r remap = L x /(L y S), and none of the results 
presented here seems to be related to this time scale. Spectral 
methods are very little dissipative by nature; numerical dissipa- 
tion is kept to very small values, as can be seen by computing 
the to tal energy budget at each time step (see lLesur & L ongaretti 
2005 for a discussion of this procedure). 

2.2. Dimensionless numbers: 

All our simulations are performed in horizontally extended 
boxes, with aspect ratio L x ; L y : L z = 4 : 4 : 1, with a non 
vanishing mean vertical field of varying strength. 

Our shearing box set up is characterized by a number of di- 
mensionless numbers. Our simulations explore the dependence 
of the turbulent transport with respect to three of them 

- The magnitude of the imposed mean vertical field measured 
by 



S 2 L 2 



P = 



(10) 



where is the Alfven speed due to this mean field. This def- 
inition mimics the usual plasma B in vertically stratified 
disk obeying the vertical hydrostatic equilibrium constraint 
c s ~ £IL Z . 
- The Reynolds number, 



Re = 



SL 2 



(11) 



comparing the nonlinear advection term to the viscous dissi- 
pation. 

- The magnetic Reynolds number, 



Rm 



?3. 



(12) 



comparing magnetic field advection to the Ohmic resistivity. 
- Alternatively to one of the two Reynolds numbers, the mag- 
netic Prandtl number 



Pm — 



Rm 
~Re~ 



(13) 



Our study focuses on the dependence of the turbulent trans- 
port on B, Re and Rm, but other combinations of the Reynolds 
numbers will be used when necessary. Note however that un- 
stratified shearing boxes are also characterized by two other 
numbers that are constant for our simulations: 



- The Mach number 
Ma 



(14) 



(vanishing for our incompressible simulations) 
- A Rossby-like number q defined by 



(15) 



which is fixed at the keplerian value (3/2). The sign of q in- 
dicates the flow cyclonicity. 
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In the following we use S~ l as the unit of time and SL : as the 
unit of velocity. One orbit corresponds to r oi b = 3nS _1 - 105 _1 ; 
our simulation lengths are to be divided by about 10 to compare 
them with simulations scaled in orbit times, as is commonly 
done by other workers in the field. For simplicity, in what fol- 
lows, we keep the same notation for dimensionless and dimen- 
sional quantities. 

Our dimensionless transport coefficient a is defined as 



a - 



(b X by - V x Vy) 

~SH 2 ' 



(16) 



where the average is taken over the volume of the simulation 
box. Within a factor of order unity, this is analogous to the more 
common prescription for the turbulent viscosity v, = ac s H if one 
identifies L, with the disc thickness H. 



channel modes that is required in order to analyze the simula- 
tions presented in the next section. 

3.1. Dispersion relation 

The physics of viscous and resis ti ve ch annel modes has been 
examined in iLesur & Long arettil (I2007I) and their pr o perties 



have been characterize d in detail in iPessah & Chad (2008) 
dSano & Mivamal[T999l have also explored the role of resistiv- 
ity on MRI in the absence of viscosity). Since non axisymmetric 
MR1 mod es are transient ly growing structures in non ideal MHD 
(Balb us & Hawlevl [l992). the axisymmetric modes give a good 
grasp of the linear stability properties of the shearing box. 



3. Channel mode physics: a summary of relevant 
information: 

As recalled in the introduction, it has been noted in a number of 
earlier MRI simulations in a shearing box with a mean vertical 
field that channel modes constitute somewhat recurrent patterns 
of the turbulent flow, and appear to be more prominent at the 
maxima in the fluctuations of the turbulent transport. Channel 
modes do transport angular momentum with roughly the right 
order of magnitude if their amplitude is comparable to the rms 
fluctuation in the field close to a maximum of the transport. 

This has suggested a picture of turbulent transport in shear- 
ing boxes with a mean vertical field that is sketched on Fig.Q] 



Channel 

Linear inst. = . mode ^ Sec. inst. 
source sink 



Small amplitude 
random motions 



Parasitic 
mode 



Nonlinear transfer 



Fig. 1. Cartoon of the process that may be at the origin of the 
saturation of the transport in shearing box simulations of MRI 
with a net vertical field (see text). 

According to this picture, the channel modes linearly grow 
from random noise; at some amplitude, their growth is halted 
by a secondary instability (parasitic mode) which destroys the 
channel mode. Presumably, the parasitic modes themselves de- 
cay into small scale turbulence, which may then produce the seed 
for the random fluctuations out of which the channel mode grows 
in the first place. In this scenario, the channel mode(s) would be 
responsible for most of the transport, in particular near maxima 
of the transport fluctuations. 

One of the objectives of this paper is to examine the rele- 
vance of this picture of MRI turbulence. Even if such a scenario 
is not generic (it depends directly on the presence of a mean ver- 
tical field), if confirmed, it might provide an interesting lead to 
analyze different situations. Another related objective it to as- 
sess the relevance of this scenario to the question of the Prandtl 
number dependence of MRI-driven transport. To this effect, we 
gather here the relevant pieces of information on the physics of 




Fig. 2. Maximal growth rate cr M /S as a function of the resistive 
dissipation, for the various values of the Prandtl number used in 
the simulations. Also shown are asymptotic analytic approxima- 
tions (see text). 

We recall here the dispersion relation of theses modes and a 
number of other features that will be of use in the next section. 

Looking for solutions of the linearized equations of motions 
in the form v = vi exp(<xf - ik y y — ikz) and b — bi exp(crf - ikz) 
leads to the following fourth order dispersion relation: 



cr 4 +2k 2 cr 3 (r] + v) + <r 2 {a + k 4 (r] 2 + v 2 + 4t]v) + bj 
r(^2k 6 (j]v 2 + v?7 2 ) + ak 2 (v + if) + Ibrjk 2 ^ 



+v 1 rj 1 k A + avrjk 4 + brfk - c = 



with 



a = 2k 2 V 2 A , 
b = k 2 j 2 , 

c = klvldSlSy 2 -klvl), 



(17) 



(18) 
(19) 
(20) 



and where k = [20(20 - S)] ' is the epicyclic frequency, Va = 
Bq is the Alfven speed based on the imposed mean vertical fielcfl 
and y 2 = k 2 /(k 2 + k 2 ). Channel modes have k v = so that y = 1, 



We measure the magnetic field in units of Alfven speed. 
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Fig. 3. The dimensionless turbulent transport a as a function of time (in units of S _1 ) for a vertical magnetic field strength /3 = 1000 
and various levels of viscous and resistive dissipation Re and Rm. Insets on a given line correspond to constant magnetic Prandtl 
number, Pm = 1/4, 1 and 4 from top to bottom. The Reynolds number Re is constant on rows (800, 1600, 3200, and 6400 
respectively. The transport averages are computed over the last 400 shear times of the simulations. The blue line represents the 
Reynolds stress, the green line the Maxwell stress, and the red line the total stress (these last two are almost indistinguishable). 



This equation is used later on to compute the channel mode 
growth rates in the conditions of our numerical simulations. It is 
most conveniently solved in dimensionless form, with S as unit 
of time and Va/S as unit of length. We introduce <x* = cr/S , 
k* = k/S , k* = VaIc/S; defining the viscous (E v ) and resistive 
(£„) Elsasser-like numbers by 



(21) 
(22) 



one has v* = l/E v and rf = 1/E n . Note that with a power- 
law velocity profile in a disk (Q oc r~ q ), Q* = 1/q and a:* = 
[2(2 - q)] l ^ 2 /q. Our definitions of the Elsasser numbers (and 
other dimension less quantities) differ fro m the usual ones (A v 
and A n , see e.g. iPessah & Chan1l2008l and lPessah & Goodman! 
2009) by a factor S/Q.. This choice has several motivations: S 
rather than Q is the relevant time-scale of linear growth rates; 
this definition is more consistent with our previous choice of 
units; and it makes the connection between the Elsasser and 
Reynolds numbers simpler. 



E v = 


V 2 

A 


Re 


Ay 




Sv ' 






E n = 


V 2 

v A 


Rm 


= A, 






~~ T 





3.2. Stability limits 

In the dissipation-free limit, the dispersion relation can be solved 
exactly. The magnetic tension stabilizes the instability for k* > 
k* = (2CT) 1/2 = (2/q) l/2 . In a box of finite vertical extent L z , 
the vertical modes wavelengths are multiple of fe m ,„ = 2jt/L z . 
Axisymmetric MRI modes are therefore stabilized when k c < 
kmin, which translates into f3 > /3 C = 2n 2 S /Q 30. The maxi- 
mum growth rate in the dissipation-free regime, cr M = 5/2, is 
achieved for a wavenumber kni — k c /2. 

Resistivity and viscosity modify the marginal stability limit 
in the {Re, Rm, (3) space. Numerical investigation indicates that 
only one of the four roots is unstable for q < 2, so that the 
marginal stability limits obtains when the cr-independent term 
of the dispersion relation is equal to zero. Also, as shown by 
iLesur & Longarettil (120071) . for p > 100 and Re > 400 (two 
conditions that are satisfied in the present investigation), the 
marginal stability limit is nearly independent of the Reynolds 
number. This leads to the following expression of the marginal 
stability limit 



Rm — — 

S 



Pc 



1/2 



2n 



P 



yl3(.P-PcY 12 



3/3 



1/2 



(23) 
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Re=3200. Rm=800, {3=100 




20 40 60 80 100 

Averaging time 

Fig. 4. Transport standard deviation as a function of the time in- 
terval r used to bin the transport data, for a given run (in this 
case, Re = 3200, Rm = 800, and p = 100). This information 
is used to quantify the error in the transport from the fit oc t~ I/ ' 2 
that is expected to hold for large enough binning time t (see text 
for details). 



where the second equality applies for keplerian flows. The last 
ap proximation holds only whe n f3 is large enough, and was given 
in iLesur & Longarettil (|2007); the first expression is new and 
sensibly more precise. 

The maximum growth rate achieved by the instability is 
modified when viscosity and resistivity becom e important, i.e. 
when either E v or E n < l. lPessah & Chad d2008) give an asymp- 
totic form of the growth rate in the dissipation regime for small 
or of order unity Pm (see their Fig. 5 and their Eq. 91); it reads 
cr M = 2<ToE n /3 where <xo = 5/2 is the dissipation-free growth 
rate. The range of E n and Pm examined in this work spans 
the transition from a dissipative to a non dissipative regime for 
the linear instability. The corresponding variation of maximum 
growth rate with both E n and E v (through the Prandtl number) 
is shown on Fig. [2] along with the asymptotic expressions just 
recalled. 



4. Numerical results 

In order to explore a significant range of parameters, most of 
our runs are performed at a standard resolution (N r , N$, N z ) = 
(256, 128, 64) in real space. All simulations have an aspect ratio 
R : (f> : z = 4- : 4 : l;as argued earlier, this aspect ratio allows 
us to capture the fastest growing channel and parasitic modes. 
Three simulations have been performed at a resolution twice as 
large (512, 256, 128) in order to reach a higher Reynolds number 
(Re = 20000), and lower Prandtl numbers (down to Pm = 1/16). 
Note that spectral codes are intrinsically more resolved than fi- 
nite difference codes with the same number of points; as a con- 
sequence, one needs to adopt resolutions larger by a factor 2 
in all directions in a finite difference code such as ZEUS, Athena 
or RAMSES to obtain the same results, a point to bear in mind 
when comparing our conclusions with thos e published in the lit- 
erature with finite difference codes (see e.g. lFromang et al.82 007, 
for an explicit comparison). 

The simulations performed in this work are in the large 
Reynolds number limit (400 < Re < 20000). The relevant 
marginal stability limit depends only on Rm and is precisely 
given by Eq. d23l . When the magnetic Reynolds number is too 
close to this limit, nonaxisymmetric perturbations damp out, and 
the flow becomes nearly identical to an axisymmetric one. We 




1/16 1/8 1/4 1 

Prandtl number 

Fig. 5. The dimensionless turbulent radial momentum transport 
a as a function of the Prandtl number Pm and the field strength 
P for various Reynolds numbers, black: p = 10 2 ; blue: p = 10 3 ; 
red: p = 10 4 ; □: Re = 400; +: Re = 800; 0: Re = 1600; o; Re = 
3200; x: Re — 6400; the green starred data points correspond to 
the three more resolved runs at Re = 20000 and p = 10 3 . Power 
law fits are also shown for each value of p. 



have checked that none of the simulations discussed below is 
affected by this issue. 

With our adopted standard resolution, we can confidently 
resolve the dissipative scales up to Reynolds numbers ^ 10 4 . 
Consequently, we have explored Reynolds numbers in the range 
400 - 6400 and magnetic Reynolds numbers in the range 400 - 
12800. The magnetic Prandtl number of the simulations is either 
1 /4, 1 or 4. A smaller value [Pm = 1/16) has also been reached 
with one of our higher resolution runs with Re = 20000. 

We have performed simulations for three different vertical 
field strengths: p = 10 2 , 10 3 , and 10 4 . This range is constrained 
by two considerations. On the low end (large field), one must 
stay away from the marginal stability limit due to the magnetic 
tension. Indeed, when the incompressible MRI works too close 
to this stability threshold, the flow is strongly biased by the pres- 
ence of channel modes (an effect that incomp ressible simulations 
are m ore prone to capture, as exemplified in lLesur & Lon^arettil 
120071) . On the large end (weak field), one is limited by the fact 
that at some point, the flow is going to be dominated by the 
zero flux transport process, whose dependence on the viscous 
and resistive dissipation might be sensibly different. If one sim- 
ply assumes that one shifts from one regime to the other when 
the efficiency of the transport processes are equal, one finds that 
the transition between the mean field and zero net field regimes 
should be in the range p ~ 10 5 - 10 4 ; this limit applies to Pm > 2 
(as the zero net flux transport vanishes for lower Prandtl number 
values), and should be somewhat dependent on Pm. As discussed 
later on, for our weakest field (J3 = 10 4 ) and smallest Prandtl 
(Pm = 4), the flow also tends to become bidimensional for low 
enough Reynolds number. 

4.1. Transport time histories and standard deviation 

Figure [3] displays the behavior of the turbulent transport for the 
runs we have performed at p = 1000. The simulations are typi- 
cally run for 500 shear times, and the transport average is based 
on the last 400. In our horizontally extended simulation boxes, 
the transport fluctuations are substantially reduced with respect 
to narrower boxes. Typically, fluctuations of a factor of ~ 2 are 
observed, whereas in r : z = 1 : 1 boxes, fluctuations of an or- 
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Fig. 6. The dimensionless transport a as a function of the magnetic Reynolds number Rm (left) and Reynolds number Re (right). 
Black: ft = 10 2 ; blue: ft = 10 3 ; red: ft = 10 4 ; □: R = 400; +: R = 800; 0: R = 1600; o: R = 3200; x: R = 6400; R refers to Re on the 
left panel, and to Rm on the right one. The green stars correspond to the three more resolved runs at Re = 20000 (see text). 



der of magnitude or more are common. This reduction provides 
more precise averages for the transport, even though the averag- 
ing time is somewhat smaller than what is used for boxes with 
a R : Z = 1 : 1 aspect ratio. This aspect ratio dependence is re- 
lated to the prominence of channel modes in the large transport 
bursts observed in the narrow boxes dBodo e t al. 2008), a feature 
related to the fact that width of the narrow box is comparable to 
the correl ation length of turbulent fluctuations in the horizontal 
direction dGuan et alJl2009t) . 

The transport time history can be used to quantify the error 
in the determination of a in the following way. From the raw 
time history, one can define a series of binned time histories, 
with a binning time r ranging from 1 to 100 shear timefl From 
these binned time histories, one can define a transport standard 
deviation in the usual way, i.e., <r a = - a) 2 /N] ' 2 where 

N is the number of bins, a, itself being the average value of the 
transport in bin ;. The resulting dependence cr a (r) is shown on 
Fig-0for one of our runs. Two regimes can be distinguished. For 
t < 10-20 S~ l , the deviation decreases sharply with t; this is 
seen directly on Fig. [3] where the transport typical variation time 
scale is precisely in this range. For large values of r, the transport 
is less and less correlated with itself. Consequently, one expects 
that it should behave more or less as a random walk, so that the 
standard deviation should scale like t -1 ^ 2 . A fit of this type is 
performed on Fig. |H and appears to represent reasonably well 
the trend, in spite of the crudeness of the approximation. 

From this approximation we find that the relative error in 
the computation of the transport over 400 shear times is in the 
range 3 - 10%. There is no clear tendency for the larger values of 
the relative error to correspond to higher transport, and no clear 
trend with either the magnetic Prandtl or magnetic Reynolds 
number. On average, cr a ja ~ 5%. 

4.2. Physical dissipation and transport 

The dependence of the transport on the dissipation for various 
field strengths is presented on Fig.[5]for all our standard resolu- 



4 For t < S 1 , the transport statistics depends very little on the bin- 
ning time, while for r > 1005 _1 , the statistics is too poor to be mean- 



ingful. 



tion runs. This figure sign ificantly extends the results reported in 
dLesur & Long aretti 2007). Several trends can be observed: 

- For a given field strength, the dependence on the magnetic 
Prandtl number Pm seems to follow an approximate power 
law. For comparison with our previous work, power-law fits 
are computed. The exponents are 0.3, 0.7 and 0.8 for the 
three values of ft, in increasing order. 

- The three more resolved runs confirm this trend: the trans- 
port is substantially larger than other runs at the same field 
strength (ft = 10 3 ), but the Pm dependence is consistent with 
the other runs at the same ft; this Pm dependence does not 
seem to saturate even though a lower Prandtl number has 
been achieved (1/16). 

- At a given Prandtl number, the transport is decreasing with 
decreasing field strength (increasing ft). This trend is ex- 
pected from the sc aling a oc /T 1/ 2 that was first pointed out in 
the initial work of lHawlev et al.l (fl9 95). However, the varia- 
tion of the index of the Prandtl number dependence with field 
strength just discussed induces deviations from this scaling. 

- There is a weaker, but systematic and significant increase of 
the transport with increasing Reynolds number at any given 
field strength and Prandtl number. This effect is real: for most 
of our runs, this increase is larger than the standard deviation 
in the transport, as quantified in the previous subsection. It is 
also larger for the smaller Prandtl number values. This indi- 
cates that the Prandtl number does not capture all the physics 
of the correlation between transport and physical dissipation; 
this point is further discussed below. 

Our previous investigation was limited to ft — 100. In the 
present work, the Prandtl number dependence of the transport 
for this field strength is consistent with our earlier findings. 
However, the transport observed here is reduced by a factor ~ 2; 
this is a direct consequence of the reduced role played by the 
channel mode in our horizontally extended simulation boxes. 

For the lowest Prandtl number (1 /4) and lowest field strength 
(ft = 10 4 ), only one point is reported in the graph. Our other 
runs for this parameter have lower Reynolds numbers, and are 
too close to the linear stability threshold of Eq. ( 1231 to sustain 
full 3D turbulent motions. These runs show 2D or quasi-2D be- 
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Fig. 7. Composite figures showing the dimensionless growth rate cr/S of the channel mode for most of the various runs performed, 
as a function of the dimensionless wavenumber Va^JS. Each subplot is labelled by its Reynolds number and field strength (J3). 
The various curves in each subplot correspond to Pm = 1/4, 1, 4, respectively (the growth rate increases with increasing Prandtl 
number). The crosses indicate the modes that are represented in our simulation boxes, as a result of the boundary conditions (see 
text for a discussion). 



havior, with very different transport efficiency and behavior. The 
data point we have retained might still be weakly affected by 
such effects. 

As pointed out above, Fig. [5] indicates that the spread with 
Reynolds number at a given Prandtl number increases with de- 
creasing Prandtl number. In fact, two different regimes can be 
noted, one for Pm < 1 and one for Pm = 4. 

At Pm = 4, the transport seems to be only weakly depen- 
dent on Re (or Rm), at least for large enough Reynolds number: 
the transport increases by 10% to 50% (depending on {$) while 
the Reynolds number is multip lied by a factor of 4. Th is trend 
can also be found in the work of Simon & Hawlev (2009), where 
the MRI turbulent transport in presence of a toroidal field is in- 
vestigated with more emphasis on the Pm > 1 regime. Their 
Figs. 6 and 7 show that, for Pm = 2 and 4 at least (the only ones 
with enough data in the Pm > 1 regime), the transport increases 
steadily with the Reynolds number for Re < 1000 and much 
more weakly for Re > 1000. 

On the contrary, the spread in Reynolds number for Pm < 1 
is substantial, and systematic. Such a spread was not detected in 
our earlier investigation, due to the larger fluctuations in trans- 
port related to the box aspect ratio, as discussed earlier. In fact, 



this dispersion seems to be an effect of the magnetic Reynolds 
number. To illustrate this point, the transport is represented on 
Fig- as a function of Rm (left panel) and Re (right panel), for 
Pm < 1; the colors describe different field strengths (fi — 10 2 to 
10 4 from top to bottom). The statistics in the number of points 
at any given Re or Rm is rather low; however, it appears quite 
clearly that the dispersion of the points at any given Reynolds 
number is substantially larger in Re (with varying Rm) than 
in Rm (with varying Re). The largest Reynolds number data 
strongly support this conclusion. Furthermore, the fits^ of the 
transport as a function of Rm indicate the Rm dependence of the 
transport for Pm < 1 is very similar to its Pm dependence as 
shown on Fig. [5] This strongly suggests that the Pm dependence 
observed on this figure is in fact mostly a Rm dependence for 
Pm < 1. Including the Pm = 4 data destroys this correlation, 
which strengthens the idea that there are two regimes, depend- 
ing on the Prandtl number (a feature that may be related to the 
existence of a transition around Pm = 2 in zero net flux shear- 



5 The Re = 20000 data points have not been included in this fit to 
make the comparison between the two dependences in the same condi- 
tions. 
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Fig. 8. The relation between dimensionless turbulent radial momentum transport a and the dimensionless growth rate cr/S . There 
is no correlation at the largest ft and Re (left) while there is one for weaker field strengths (right). In the left panel, the transport and 
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Pm =1/4 run have been used to normalize the data. Black: ft 
= A. The three resolved runs are also shown (green stars). 
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ing bo x simulations). The relevant results of I Simon & Hawlevl 
(2009); although less detailed, are consistent with these findings 
(see their Fig. 7). 

5. Role of channel and parasitic modes 

5.1. Linear physics and turbulent transport 

iLesur & Longaretti d2007l) concluded that there was no direct 
connection between the Prandtl dependence of MRI-driven tur- 
bulent transport and the linear growth rate of channel modes. 
However, the transport was substantially less precisely deter- 
mined than in the present investigation, and the issue is worth 
reinvestigating. 

To this effect, one needs to compute the linear instability 
growth rates for the various runs we have performed. We fo- 
cus on the growth rate of the channel mode; for one thing, it is 
a guideline for the stability behavior of all unstable modes, and 
for another, there is a particular interest in the behavior of this 
specific mode. The resulting growth rates are numerically com- 
puted from the dispersion relation recalled in section I3.ll and 
represented in Fig. [7] The important features to note from this 
figure are the following: 

- The largest growth rate in the continuum limit is always rep- 
resented or very nearly represented in the simulation, in spite 
of the discreteness imposed by the vertical boundary condi- 
tion (the corresponding modes are shown by crosses on the 
figure). This ensures that the fastest growing channel mode 
is always present in our simulations. 

- For the largest field strength (ft = 100), the growth rates 
of all our simulations is very close to the ideal MHD limit 
(<x = S 12). This is especially true of the two largest Reynolds 
numbers (Re = 3200 and Re = 6400). Note in particular for 
these runs that the growth rate changes by at most 10% while 
the transport varies by a factor of more than 3 (see left panel 
of Fig.. [8). As a consequence the existence of the Prandtl and 
magnetic Reynolds number dependences of the transport is 
not related to the variation of the growth rate with physical 
dissipation. 

- For ft = 1000, and for the two largest Prandtl numbers 
(Pm = 1 and Pm = 4) and largest Reynolds numbers 
(Re = 3200 and Re = 6400), the growth rate is again almost 



independent of the dissipation. As there is only one captured 
channel mode for the /3 = 100 vs four or five for ft = 1000, 
the Pm, Rm dependence is not related to the number of cap- 
tured channel modes in the simulation. 

There is a sharp behavior difference of the transport with re- 
spect to the linear MRI growth rate between the nearly ideal sim- 
ulations performed at ft — 100 and the lower ft values, as is ap- 
parent when comparing the left and right panels of Fig. [8] In fact, 
the whole spectrum of possible behaviors is covered here, from 
substantial variation of the transport with very weakly varying 
growth rates (left panel), or conversely substantial variation of 
the growth rate while the transport is only weakly or mildly 
changing (Pm = A data on the right panel) and intermediate situ- 
ations where the two are more strongly correlated (the other data 
points on the right panel). Ther efore, although this situat ion is 
more contrasted than indicated in lLesur & Longarettil(l2007l) . the 
conclusion stands: there is no clear and direct relation between 
the effects of dissipation on linear physics and the nonlinear tur- 
bulent saturation mechanism. 

5.2. Quasi-linear physics and turbulent transport: 

As indicated in the introduction and in section [3j it is often as- 
sumed that the channel mode is responsible for most of the MRI- 
driven turbulent transpor t in the presence of a net ve rtical mag- 
netic field. In particular Pessah & Goodman (2009) argue that 
the saturation of the channel mode by its parasites may help us 
to understand the transport. In this section, we make use of our 
numerical results to investigate these issues. 

We first disprove the first statement: the channel mode does 
not clearly dominate transport, especially at the maxima of the 
transport, as should in the scheme depicted on Fig. [1] This is 
shown explicitly in Fig. [9] A one-dimensional Fourier trans- 
form of both the Maxwell and Reynolds stress along the az- 
imuthal direction is shown in this figure. These transforms are 
performed at a typical minimum and maximum of transport in 
the transport history of two of our runs, namely the runs with 
Re = Rm = 6400 at ft = 100 and ft = 10 4 , respectively. The 
transform is shown in terms of the wavenumber n = k^L^/ln; 
n — corresponds to the transport due to the axisymmetric 
modes (channel modes) present in the simulations, for the vari- 
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transport maximum in the time history, while the lower ones describe a typical minimum. Both the Reynolds and Maxwell stress 
are shown. The transport due to non-axisymmetric modes always dominates the transport due to the axisymmetric ones (channel 
modes). See text. 



ous captured k z . The other wavenumbers (« + 0) correspond to 
non-axisymmetric modes. Channel modes dominate over other 
individual wavenumbers for ft = 100. However, the cumulated 
contribution of non-axisymmetric modes is dominant by a factor- 
s' at the maximum of transport for this ft, and does not exceed 
the non-axisymmetric contribution by more than ~ 30% at min- 
imum transport. For ft = 10 4 , the cumulated non-axisymmetric 
are always dominant by a factor ~ 6. For ft = 10 3 (not shown), 
the situation is somewhat intermediate between the two extreme 
values presented in the figure, but the non-axisymmetric contri- 
bution is again always dominant by a factor several. 

In summary, except for minimum transport for the largest 
field strength where axisymmetric and non-axisymmetric con- 
tributions to the transport are comparable, the non-axisymmetric 
contribution is always larger by a factor ~ 5. As a result, it seems 
necessary at least to include non-axisymmetric modes in the pro- 
cess sketched on Fig.Q] 

Furthermore, the analysis of Pess ah & Goodman! (l2009h is 
inconsistent with our results, both from a qualitative and quan- 
titative point of view. Our simulations always capture the fastest 
growing channel mode, or very nearly so (see Fig. [7]). Due to the 
choice of aspect ratio, the fastest growing parasitic modes are 
always captured as well. Finally, for ft = 100 and Re = 3200 
and 6400 (even 1600 for the two larger Pm values), the vari- 
ation of the transport wit h Pm is quite substantial, whereas 
Pessah & Goodma U (120091) conclude that in these conditions, 
the transport is independent of dissipation (see the left plot in 
Fig. 2 of their paper; the three Reynolds numbers values men- 
tioned correspond to the three largest Elsasser numbers of their 
graph). 

6. Discussion and conclusion 

Perhaps the most significant new result of this work, disclosed 
on Figs. [5] and [6] is the existence of a double regime separated 



by a critical magnetic Prandtl number Pm c ~ 1 . For Pm < Pm c , 
at a given field strength, the transport correlates mostly with Rm; 
for Pm > Pm c , the transport seems to depend mostly on Pm 
and only weakly on either Re or Rm (once Re > 10 3 ), although 
a larger number of Pm values need to be probed on this issue. 
It is tempting to assume Pm c =t 2, as this is the critical value 
for the zero mean field problem, but this identification requires 
further work to be substantiated. The identification of this double 
dissipation regime was made possible by the increased accuracy, 
with respect to our previous work, in the determination of the 
transport averages. In the small Prandtl regime, in contrast to 
the large one, our most resolved simulations show no sign of 
convergence with respect to dissipation, although values of Rm 
up to 20000 have been reached. 

The role of linear physics and parasitic modes on transport 
properties has also been investigated, and the major results on 
these questions can be summarized as follows: 

- The existence of the dependence of turbulent transport on 
dissipation is not related to the role of dissipation on linear 
modes growth rates. 

- The existence of the correlation of transport with dissipation 
is not related to the number of channel modes captured in the 
simulations. 

- Saturation of the channel modes growth by parasitic modes 
is not the responsible for correlation of transport with dis- 
sipation. This is hardly surprising as both the channel and 
parasitic modes are large scale modes, that are expected to 
be little affected by dissipation, whereas the trends of the 
transport with dissipation are substantial. 

- Furthermore, the transport is usually not dominated by ax- 
isymmetric modes; these modes contribute at most at the 
same level as non-axisymmetric modes for the strongest 
fields investigated, and have negligible contributions for the 
weakest fields. 
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Note that all these results were obtained while the fastest 
growing channel and parasitic modes were always captured, so 
that they do not depend on limitations on this front. 

It is worth pointing out that our results do not totally dis- 
qualify a saturation of the unstable linear modes by the parasitic 
modes; only the relevance of this process to the relation between 
transport and physical dissipation has been disproved. 

The role of the Prandtl number disclosed in this investiga- 
tion makes more physical sense than a direct dependence of the 
transport on Pm. Indeed, it is very likely that the critical value 
Pm c relates to the switching of the magnetic and kinetic dissi- 
pations scales k n and k v in Fourier space: for Pm < Pm c (resp. 
Pm > Pm c ), k n < k v (resp. k, ; > k v ). For Pm «: 1, k n <sc k v , the 
flow at scales k > or » k n in Fourier space is therefore purely 
hydrodynamic. As information flows from small k to large k in 
this scale range, the flow should be independent of Re, so that in 
this regime a = a(J3, Rm) is expected. The transport may even 
become independent of Rm at large enough Rm in the absence 
of backreaction of the small scales on the large ones, and if en- 
ergy exchanges are not strongly nonlocal in Fourier space. We 
merely note here that in the context of homogeneous, isotropic, 
incompressible MHD turbulence, all energy transfers are direct 
(from small k to large k); exchanges betwe en magnetic and ki - 
netic energy are nonlocal in Fourier space dAlexakis et alj2007l) . 
altho ugh this nonlocality seems to be of finite (albeit large) ex- 
tent dAluie & Evink 2009). The situation is more complex in the 
large Prandtl regime, as non-local transfers and small scale dy- 
namo action might take place at scales smaller than the viscous 
dissipation scale, and further investigations are required to char- 
acterize the properties of the small dissipation limit, in particular 
concerning the locality of transfers. 

The behavior of MRI-driven turbulence with respect to these 
various issues will be reported elsewhere, through the analysis of 
energy transfers in Fourier space. More generally, investigating 
the possible independence of transport with respect to dissipa- 
tion in the Pm <s Pm c or Pm » Pm c limits requires to achieve 
a double scale separation: one must first separate in Fourier 
space the smallest injection scales (the scales which contribute to 
the Reynolds and Maxwell stresses) from the largest dissipation 
scales, and then the two dissipation scales themselves. This is 
extremely demanding in terms of resolution. Our simulations in 
the small Prandtl regime are consistent with a separation of dis- 
sipation scales nearly achieved. As a first step at small Prandtl 
(the most critical regime for YSO disks), one can achieve one 
or the other of the two scales separation. Such investigations are 
underway and will be reported elsewhere. 

These questions are critical for astrophysical disk dynamics, 
where the Reynolds numbers are always very large (Re > 10 10 , 
while the Prandtl number spans very small (< 10~ 5 , YSO disks) 
to very large values (> 10 5 , AGN disks). Addressing these is- 
sues requires to overcome a number of limitations in the sim- 
ulations. Besides the questions of scale separation and nonlo- 
cality of transfers in Fourier space just mentioned, the role of 
compressibility and vertical stratification on these results must 
be quantified. These questions can certainly be investigated in 
shearing boxes, but global simulations are probably still largely 
out of reach, due to the resolution demand. Also (and probably 
as much importantly), different dissipation regimes must be an- 
alyzed in the context of YSOs, namely the ambipolar and Hall 
regime, which are known to have substantial effect on the linear 
development of the MRJ, while being relevant for large fractions 
of YSO disks. 
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